Source code and data found at https://github.com/maRce10/budgie_inbre_stress_vocal_output

Load packages

## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2",
    "kableExtra", "knitr", "formatR", "MASS", "sp", "GGally", "brms",
    "lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes",
    "cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace")

source("~/Dropbox/R_package_testing/sketchy/R/load_packages.R")
load_packages(x)

source("~/Dropbox/R_package_testing/brmsish/R/html_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
source("~/Dropbox/R_package_testing/brmsish/R/read_summary.R")

Functions and global parameters

opts_knit$set(root.dir = "..")

# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE,
    tidy = TRUE)

read_excel_df <- function(...) data.frame(read_excel(...))


# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")


standard_error <- function(x) sd(x)/sqrt(length(x))

cols <- viridis(10, alpha = 0.7)


# print results brms models
summary_brm_model <- function(x, gsub.pattern = NULL, gsub.replacement = NULL,
    xlab = "Effect size") {

    summ <- summary(x)$fixed
    fit <- x$fit
    betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE)
    hdis <- t(sapply(betas, function(y) hdi(fit@sim$samples[[1]][[y]])))
    summ_table <- data.frame(summ, hdis, iterations = attr(fit, "stan_args")[[1]]$iter,
        chains = length(attr(fit, "stan_args")))
    summ_table <- summ_table[rownames(summ_table) != "Intercept",
        c("Estimate", "Rhat", "Bulk_ESS", "l.95..CI", "u.95..CI",
            "iterations", "chains")]

    summ_table <- as.data.frame(summ_table)
    summ_table$Rhat <- round(summ_table$Rhat, digits = 3)
    summ_table$CI_low <- round(unlist(summ_table$l.95..CI), digits = 3)
    summ_table$CI_high <- round(unlist(summ_table$u.95..CI), digits = 3)
    summ_table$l.95..CI <- summ_table$u.95..CI <- NULL

    out <- lapply(betas, function(y) data.frame(variable = y, value = sort(fit@sim$samples[[1]][[y]],
        decreasing = FALSE)))

    posteriors <- do.call(rbind, out)
    posteriors <- posteriors[posteriors$variable != "b_Intercept",
        ]

    if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
        posteriors$variable <- gsub(pattern = gsub.pattern, replacement = gsub.replacement,
            posteriors$variable)

    gg <- ggplot(data = posteriors, aes(y = variable, x = value, fill = stat(quantile))) +
        geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = HDInterval::hdi,
            vline_linetype = 2) + theme_classic() + labs(x = xlab,
        y = "Predictor") + xlim(range(c(min(summ_table[, "CI_low"]),
        max(summ_table[, "CI_high"])), 0)) + geom_vline(xintercept = 0,
        col = "red") + scale_fill_manual(values = c("transparent",
        "lightblue", "transparent"), guide = "none") + ggtitle(x$formula)

    if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
        rownames(summ_table) <- gsub(pattern = gsub.pattern, replacement = gsub.replacement,
            rownames(summ_table))

    summ_table$Rhat <- ifelse(summ_table$Rhat > 1.05, cell_spec(summ_table$Rhat,
        "html", color = "white", background = "red", bold = TRUE,
        font_size = 12), cell_spec(summ_table$Rhat, "html"))

    signif <- summ_table[, "CI_low"] * summ_table[, "CI_high"] > 0

    df1 <- kbl(summ_table, row.names = TRUE, escape = FALSE, format = "html",
        digits = 3)

    df1 <- row_spec(kable_input = df1, row = which(summ_table$CI_low *
        summ_table$CI_high > 0), background = adjustcolor(cols[9],
        alpha.f = 0.3))

    df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover",
        "condensed", "responsive"), full_width = FALSE, font_size = 12)

    print(df1)

    print(gg)

}
# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")

# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd",
    "freq.median", "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent",
    "time.ent", "entropy", "meandom", "mindom", "maxdom", "dfrange",
    "modindx", "meanpeakf")]

sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])

sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year,
    sep = "-")

sels$date <- as.Date(sels$date, format = "%d-%b-%Y")

# acoustic measurements
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")

indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")

indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")

# group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID),
# function(x) { X <- group_ovlp[group_ovlp$ID == x, ]
# X$overlap.to.group <- X$overlap.to.group -
# X$overlap.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] X$distance.to.group <-
# X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))

group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

Counts per individual

df <- as.data.frame(table(sels$ID))

names(df) <- c("ID", "Sample_size")

df <- df[order(df$Sample_size, decreasing = FALSE), ]

kb <- kable(df, row.names = FALSE)

kb <- kable_styling(kb, bootstrap_options = c("striped", "hover",
    "condensed", "responsive"))

print(kb)
ID Sample_size
132YGMM 6
125YGMM 12
160YGHM 16
310YGHM 21
393YGHM 25
303YGHM 37
398WBHM 41
365YGHM 44
399YGLM 46
300YGHM 47
270YGHM 64
407YGHM 97
386WBHM 100
377WWLM 108
367WWNM 119
371YYLM 149
397YGHM 175
378YYLM 195
404WBHM 196
258YGHM 223
389WWLM 230
262YGHM 306
400YGHM 360
388YGLM 377
327YYLM 404
396YBHM 455
403WBHM 566
356WBHM 761
361YGHM 767
402YGHM 849
395WBHM 854
360YGHM 900
390WBHM 935
405YBHM 975
385YBHM 1256
394YBHM 1693

Add metadata

metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")

# head(metadat)

sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"

# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID,
# sels$ID)

sels$treatment <- sapply(1:nrow(sels), function(x) {

    metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]

})


sels$treatment.room <- sapply(1:nrow(sels), function(x) {

    metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]

})


sels$round <- sapply(1:nrow(sels), function(x) {

    metadat$Round[metadat$Bird.ID == sels$ID[x]][1]

})

sels$source.room <- sapply(1:nrow(sels), function(x) {

    metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]

})

sels$record.group <- sapply(1:nrow(sels), function(x) {

    metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]

})


# add week
out <- lapply(unique(sels$round), function(x) {

    Y <- sels[sels$round == x, ]

    min_date <- min(Y$date)
    week_limits <- min_date + seq(0, 100, by = 7)

    Y$week <- NA
    for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i -
        1] & Y$date < week_limits[i]] <- i - 1

    return(Y)
})

sels <- do.call(rbind, out)


sels$cort.baseline <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


sels$cort.stress <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})

sels$stress.response <- sels$cort.stress - sels$cort.baseline

sels$weight <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


sels$breath.count <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


# remove week 5
sels <- sels[sels$week != 5, ]
agg_dat <- aggregate(selec ~ ID + week, data = sels, length)

# compare to week 1 agg_dat$call.count <-
# sapply(1:nrow(agg_dat), function(x) { baseline <-
# agg_dat$selec[agg_dat$week == 1 & agg_dat$ID == agg_dat$ID[x]]
# if (length(baseline) > 0) change <- agg_dat$selec[x] -
# baseline else change <- agg_dat$selec[x] return(change) } )

# without comparing to week 1
agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) agg_dat$selec[x])

agg_dat$selec <- NULL

agg_dat$baseline.CORT <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$cort.baseline[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$cort.baseline[sels$week == agg_dat$week[x] & sels$ID ==
        agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$stress.response <- sapply(1:nrow(agg_dat), function(x) {

    # baseline <- sels$stress.response[sels$week == 1 & sels$ID
    # == agg_dat$ID[x]]
    current <- sels$stress.response[sels$week == agg_dat$week[x] &
        sels$ID == agg_dat$ID[x]]

    # if (length(baseline) > 0 & length(current) > 0) change <-
    # mean(current) - mean(baseline) else change <- NA

    return(mean(current))
})


agg_dat$stress.CORT <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$cort.stress[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$cort.stress[sels$week == agg_dat$week[x] & sels$ID ==
        agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$weight <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$weight[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$weight[sels$week == agg_dat$week[x] & sels$ID ==
        agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$breath.rate <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$breath.count[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$breath.count[sels$week == agg_dat$week[x] & sels$ID ==
        agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$acoustic.area <- sapply(1:nrow(agg_dat), function(x) {

    area <- areas_by_week$raref.area[areas_by_week$ID == agg_dat$ID[x] &
        areas_by_week$week == agg_dat$week[x]]

    if (length(area) < 1)
        area <- NA

    return(area)
})

agg_dat$acoustic.distance <- sapply(1:nrow(agg_dat), function(x) {

    distance <- indiv_ovlp$distance.to.first.week[indiv_ovlp$ID ==
        agg_dat$ID[x] & indiv_ovlp$week == agg_dat$week[x]]

    if (length(distance) < 1)
        distance <- NA

    return(distance)
})

agg_dat$acoustic.overlap <- sapply(1:nrow(agg_dat), function(x) {

    overlap <- indiv_ovlp$overlap.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
        indiv_ovlp$week == agg_dat$week[x]]

    if (length(overlap) < 1)
        overlap <- NA

    return(overlap)
})

agg_dat$acoustic.overlap.to.group <- sapply(1:nrow(agg_dat), function(x) {

    overlap <- group_ovlp$overlap.to.group[group_ovlp$ID == agg_dat$ID[x] &
        group_ovlp$week == agg_dat$week[x]]

    if (length(overlap) < 1)
        overlap <- NA

    return(overlap)
})


agg_dat$treatment <- sapply(1:nrow(agg_dat), function(x) unique(sels$treatment[sels$ID ==
    agg_dat$ID[x]]))

agg_dat$round <- sapply(1:nrow(agg_dat), function(x) unique(sels$round[sels$ID ==
    agg_dat$ID[x]]))

Physiological parameters

breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count",
    "D14.Breath.Count", "D21.Breath.Count", "D28.Breath.Count")])

weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.",
    "D14.Bird.Weight..g.", "D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])

cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress",
    "D14.CORT.Stress", "D21.CORT.Stress", "D28.CORT.Stress")])

cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline",
    "D14.CORT.Baseline", "D21.CORT.Baseline", "D28.CORT.Baseline")])


stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round",
    "Treatment.Room")], week = breath.count$ind, breath.count = breath.count$values,
    weight = weight$values, cort.stress = cort.stress$values, cort.baseline = cort.baseline$values,
    stress.response = cort.stress$values - cort.baseline$values)

# head(stress)

stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."),
    "[[", 1), levels = c("D3", "D7", "D14", "D21", "D28"))

stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)

stress$treatment <- factor(stress$Treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

# remove 5th week
stress <- stress[stress$week != "D28", ]


stress_l <- lapply(stress$Bird.ID, function(x) {
    X <- stress[stress$Bird.ID == x, ]

    X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
    X$weight <- X$weight - X$weight[X$week == "D3"]
    X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
    X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week ==
        "D3"]
    X$stress.response <- X$stress.response  #- X$stress.response[X$week == 'D3']

    return(X)
})

stress <- do.call(rbind, stress_l)

agg_stress <- aggregate(cbind(breath.count, weight, stress.response,
    cort.baseline) ~ week + treatment, stress, mean)
agg_stress_se <- aggregate(cbind(breath.count, weight, stress.response,
    cort.baseline) ~ week + treatment, stress, standard_error)

names(agg_stress_se) <- paste(names(agg_stress_se), ".se", sep = "")

agg_stress <- cbind(agg_stress, agg_stress_se[, 3:6])

agg_stress$Week <- 1:4

bs <- 22

gg_breath.count <- ggplot(data = agg_stress, aes(x = Week, y = breath.count,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = breath.count -
    breath.count.se, ymax = breath.count + breath.count.se), width = 0.1) +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
    ncol = 3, scale = "fixed") + labs(y = "Mean change in breath\nrate (breaths/min)",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")

gg_weight <- ggplot(data = agg_stress, aes(x = Week, y = weight, fill = treatment)) +
    geom_bar(stat = "identity") + geom_errorbar(aes(ymin = weight -
    weight.se, ymax = weight + weight.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") +
    labs(y = "Mean change in \nweight (grams)", x = "Week") + theme_classic(base_size = bs) +
    theme(legend.position = "none")


gg_cort.baseline <- ggplot(data = agg_stress, aes(x = Week, y = cort.baseline,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = cort.baseline -
    cort.baseline.se, ymax = cort.baseline + cort.baseline.se), width = 0.1) +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
    ncol = 3, scale = "fixed") + labs(y = "Mean change in\nbaseline CORT (ng/mL)",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")


gg_stress.response <- ggplot(data = agg_stress, aes(x = Week, y = stress.response,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = stress.response -
    stress.response.se, ymax = stress.response + stress.response.se),
    width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Stress response \nCORT (ng/mL)",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")

cowplot::plot_grid(gg_breath.count, gg_weight, gg_cort.baseline, gg_stress.response)

# cowplot::ggsave2(filename =
# './output/physiological_parameters_tru_time_70dpi.jpeg')
# cowplot::ggsave2(filename =
# './output/physiological_parameters_tru_time_300dpi.jpeg', dpi
# = 300)

Stats

Models: Predicted physio measure ~ treatment + week (continuous) + IndRandom

Variables (Difference from Week 1): weight, BR, baseline CORT, Stress CORT, Stress Response

responses <- c("baseline.CORT", "stress.response", "stress.CORT",
    "weight", "breath.rate")

predictors <- c("~ treatment + week + (1|ID) + (1|round)")

formulas <- expand.grid(responses = responses, predictors = predictors,
    stringsAsFactors = FALSE)

vars_to_scale <- c(responses, "week")

# remove week 1
sub_agg_dat <- agg_dat[agg_dat$week != 1, ]

for (i in vars_to_scale) sub_agg_dat[, vars_to_scale] <- scale(sub_agg_dat[,
    vars_to_scale])

physio_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- sub_agg_dat[!is.na(sub_agg_dat[names(sub_agg_dat) ==
        formulas$responses[x]]), ]
    sub_dat

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 20000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 4), silent = TRUE)
    return(mod)
})


names(physio_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(physio_models, "./data/processed/physiological_response_models.RDS")
physio_models <- readRDS("./data/processed/physiological_response_models.RDS")

ggl <- list()
for (x in 1:length(physio_models)) {

    ggl[[length(ggl) + 1]] <- html_summary(physio_models[[x]], gsub.pattern = "b_treatment|b_",
        gsub.replacement = "", highlight = FALSE, remove.intercepts = TRUE,
        model.name = names(physio_models)[x])
    print(ggl[[length(ggl)]])
}

baseline.CORT ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 Intercept-student_t(3, 0.1, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) baseline.CORT ~ treatment + week + (1 | ID) + (1 | round) 20000 4 1 10000 852 0 4603.542 2994.716 905461639
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress 0.885 0.125 1.656 1.001 5810.772 11998.851
MediumStress 0.195 -0.632 1.009 1.001 4603.542 2994.716
week -0.012 -0.137 0.111 1 18860.740 19591.347

stress.response ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 Intercept-student_t(3, -0.1, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) stress.response ~ treatment + week + (1 | ID) + (1 | round) 20000 4 1 10000 291 0 17259.93 7182.645 1959704000
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress -0.155 -0.742 0.443 1 21183.71 24652.050
MediumStress -0.078 -0.674 0.537 1 17259.93 7182.645
week -0.108 -0.290 0.075 1 42051.84 25878.502

stress.CORT ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 Intercept-student_t(3, 0.1, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) stress.CORT ~ treatment + week + (1 | ID) + (1 | round) 20000 4 1 10000 1382 0 1215.206 295.917 1998039246
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress -0.252 -1.032 0.534 1 6526.781 14743.384
MediumStress 0.095 -0.790 0.949 1.003 4011.007 3739.939
week -0.080 -0.195 0.031 1.003 1215.206 295.917

weight ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 Intercept-student_t(3, 0.1, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) weight ~ treatment + week + (1 | ID) + (1 | round) 20000 4 1 10000 660 0 7094.293 12124.98 815174156
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress -0.413 -1.232 0.405 1.001 7094.293 13001.57
MediumStress -0.190 -1.090 0.710 1.001 7498.298 12124.98
week -0.086 -0.211 0.041 1.001 11368.583 23277.30

breath.rate ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 Intercept-student_t(3, -0.4, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) breath.rate ~ treatment + week + (1 | ID) + (1 | round) 20000 4 1 10000 536 0 17035.19 9553.508 1346894033
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress 0.219 -0.171 0.617 1 17035.19 11812.759
MediumStress 0.110 -0.317 0.544 1 18478.62 9553.508
week -0.771 -0.916 -0.628 1 35543.50 23097.254

names(ggl) <- names(physio_models)

 

Barplot and effect sizes side-by-side

cowplot::plot_grid(gg_breath.count + theme_classic(base_size = 10) +
    theme(legend.position = "none"), ggl[[grep("breath", names(ggl))]])

cowplot::plot_grid(gg_weight + theme_classic(base_size = 10) + theme(legend.position = "none"),
    ggl[[grep("weight", names(ggl))]])

cowplot::plot_grid(gg_cort.baseline + theme_classic(base_size = 10) +
    theme(legend.position = "none"), ggl[[grep("base", names(ggl))]])

cowplot::plot_grid(gg_stress.response + theme_classic(base_size = 10) +
    theme(legend.position = "none"), ggl[[grep("response", names(ggl))]])

Takeaways

  • Breath rate decreases gradually with time across after the first week

  • Stress response is higher in “high stress” birds compared to first week

 

Acoustic space projection

t-SNE

scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median",
    "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent",
    "entropy", "meandom", "mindom", "maxdom", "dfrange", "modindx",
    "meanpeakf")])

tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE,
    max_iter = 5000)

saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")

Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")

sels <- data.frame(sels, Y)

sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) +
    geom_point() + labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) +
    theme_classic(base_size = 25) + guides(colour = guide_legend(override.aes = list(size = 10)))

Behavioral parameters

agg_call.count <- aggregate(cbind(call.count, acoustic.overlap.to.group) ~
    week + treatment, agg_dat, mean)

agg_behav <- aggregate(cbind(acoustic.area, acoustic.overlap) ~ week +
    treatment, agg_dat, mean)

agg_call.count_se <- aggregate(cbind(call.count, acoustic.overlap.to.group) ~
    week + treatment, agg_dat, standard_error)

agg_behav_se <- aggregate(cbind(acoustic.area, acoustic.overlap) ~
    week + treatment, agg_dat, standard_error)

agg_behav_se <- merge(agg_call.count_se, agg_behav_se, all = TRUE)

names(agg_behav_se) <- paste(names(agg_behav_se), ".se", sep = "")

agg_behav <- merge(agg_call.count, agg_behav, all = TRUE)

agg_behav <- cbind(agg_behav, agg_behav_se[, 3:6])

bs <- 22

agg_behav$treatment <- factor(agg_behav$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

gg_call.count <- ggplot(data = agg_behav, aes(x = week, y = call.count,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = call.count -
    call.count.se, ymax = call.count + call.count.se), width = 0.1) +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
    ncol = 3, scale = "fixed") + labs(y = "Vocal output", x = "Week") +
    theme_classic(base_size = bs) + theme(legend.position = "none")

gg_acoustic.area <- ggplot(data = agg_behav, aes(x = week, y = acoustic.area,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.area -
    acoustic.area.se, ymax = acoustic.area + acoustic.area.se), width = 0.1) +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
    ncol = 3, scale = "fixed") + labs(y = "Change in vocal diversity",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")

gg_acoustic.overlap <- ggplot(data = agg_behav, aes(x = week, y = acoustic.overlap,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.overlap -
    acoustic.overlap.se, ymax = acoustic.overlap + acoustic.overlap.se),
    width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Vocal Plasticity",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")

gg_acoustic.overlap.group <- ggplot(data = agg_behav, aes(x = week,
    y = acoustic.overlap.to.group, fill = treatment)) + geom_bar(stat = "identity") +
    geom_errorbar(aes(ymin = acoustic.overlap.to.group - acoustic.overlap.to.group.se,
        ymax = acoustic.overlap.to.group + acoustic.overlap.to.group.se),
        width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Vocal convergence",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")

cowplot::plot_grid(gg_call.count, gg_acoustic.area, gg_acoustic.overlap,
    gg_acoustic.overlap.group)

# cowplot::ggsave2(filename =
# './output/vocal_parameters_tru_time_70dpi.jpeg')
# cowplot::ggsave2(filename =
# './output/vocal_parameters_tru_time_300dpi.jpeg', dpi = 300)

Stats

Model: Predicted behavior ~ treatment + week (continuous) + IndRandom

Variables: # calls, Distance moved from self in first week, Overlap to original acoustic space, Match to group repertoire, Maybe overall size of acoustic space

Do as comparison to week one using rarefacted calls and kernel density

responses <- c("call.count", "acoustic.area", "acoustic.overlap",
    "acoustic.overlap.to.group")

predictors <- c("~ treatment + week + (1|ID) + (1|round)")

formulas <- expand.grid(responses = responses, predictors = predictors,
    stringsAsFactors = FALSE)

vars_to_scale <- c(responses, "week")


for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[,
    vars_to_scale])

behav_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
        ]

    # remove week 1
    if (!grepl("count|group", formulas$responses[x]))
        sub_dat <- sub_dat[sub_dat$week != 1, ]

    mod <- brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 50000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9,
            max_treedepth = 15), chains = 4, prior = c(prior(normal(0,
            5), "b"), prior(normal(0, 10), "Intercept"), prior(student_t(3,
            0, 10), "sd"), prior(student_t(3, 0, 10), "sigma")))

    return(mod)
})


names(behav_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(behav_models, "./data/processed/behavioral_response_models.RDS")
behav_models <- readRDS("./data/processed/behavioral_response_models.RDS")

ggl <- list()
for (x in 1:length(behav_models)) {

    ggl[[length(ggl) + 1]] <- html_summary(behav_models[[x]], gsub.pattern = "b_treatment|b_",
        gsub.replacement = "", highlight = FALSE, remove.intercepts = TRUE,
        model.name = names(behav_models)[x])
    print(ggl[[length(ggl)]])
}

call.count ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 20) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) call.count ~ treatment + week + (1 | ID) + (1 | round) 50000 4 1 25000 1540 0 25959.37 35763.57 1897146692
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress 0.787 0.217 1.354 1 28006.98 44077.55
MediumStress 0.396 -0.183 0.967 1 25959.37 35763.57
week -0.058 -0.202 0.085 1 66399.00 56720.49

acoustic.area ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 20) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) acoustic.area ~ treatment + week + (1 | ID) + (1 | round) 50000 4 1 25000 3780 0 6759.291 5875.617 76261279
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress -0.218 -1.058 0.633 1.001 6759.291 5875.617
MediumStress -0.707 -1.558 0.151 1 9861.602 20907.370
week 0.035 -0.083 0.151 1 7041.084 9145.749

acoustic.overlap ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 20) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) acoustic.overlap ~ treatment + week + (1 | ID) + (1 | round) 50000 4 1 25000 2656 0 21041.63 31360.35 2038732099
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress 0.900 0.050 1.741 1 22453.07 37240.04
MediumStress 0.236 -0.692 1.155 1 21041.63 31360.35
week -0.169 -0.307 -0.031 1 57619.53 58721.25

acoustic.overlap.to.group ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 20) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) acoustic.overlap.to.group ~ treatment + week + (1 | ID) + (1 | round) 50000 4 1 25000 4289 0 707.655 171.241 450263308
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress 0.345 -0.303 1.016 1.001 11286.051 23797.575
MediumStress 0.361 -0.290 1.038 1.001 21649.351 39958.132
week -0.220 -0.386 -0.032 1.006 707.655 171.241

names(ggl) <- names(behav_models)

 

Barplot and effect sizes side-by-side

col_pointrange <- viridis(10)[2]

gg_coeffs <- lapply(behav_models, function(x) {

    vars <- grep("b_", posterior::variables(x), value = TRUE)
    draws <- posterior::as_draws_array(x, variable = vars)

    coef_table <- draw_summary(draws, variables = vars, probs = c(0.025,
        0.975), robust = TRUE)

    coef_table$predictor <- rownames(coef_table)
    coef_table$predictor <- gsub("b_treatment|b_", "", coef_table$predictor)
    coef_table$predictor <- gsub("Stress", " stress", coef_table$predictor)
    coef_table$predictor <- gsub("week", "Week", coef_table$predictor)
    coef_table <- coef_table[coef_table$predictor != "Intercept",
        ]

    gg_coef <- ggplot2::ggplot(data = coef_table, aes(x = Estimate,
        y = predictor)) + geom_vline(xintercept = 0, lty = 2) + ggplot2::geom_point(size = 4,
        col = col_pointrange) + ggplot2::geom_errorbar(ggplot2::aes(xmin = `l-95% CI`,
        xmax = `u-95% CI`), width = 0, col = col_pointrange) + ggplot2::theme_classic(base_size = bs) +
        ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"),
            plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "none",
            strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_blank()) +
        ggplot2::labs(x = "Effect size", y = "")

    return(gg_coef)
})

# cowplot::plot_grid(gg_call.count + theme_classic(base_size =
# 20) + theme(legend.position = 'none'),
# gg_coeffs[[grep('count', names(gg_coeffs))]], nrow = 2)
# cowplot::plot_grid(gg_acoustic.area + theme_classic(base_size
# = 20) + theme(legend.position = 'none'),
# gg_coeffs[[grep('area', names(gg_coeffs))]])
# cowplot::plot_grid(gg_acoustic.overlap +
# theme_classic(base_size = 20) + theme(legend.position =
# 'none'), gg_coeffs[[grep('overlap ~', names(gg_coeffs))]])
# cowplot::plot_grid(gg_acoustic.overlap.group +
# theme_classic(base_size = 10) + theme(legend.position =
# 'none'), gg_coeffs[[grep('group', names(gg_coeffs))]])
bs <- 10
cowplot::plot_grid(gg_call.count + theme_classic(base_size = bs) +
    theme(legend.position = "none"), gg_acoustic.area + theme_classic(base_size = bs) +
    theme(legend.position = "none"), gg_acoustic.overlap + theme_classic(base_size = bs) +
    theme(legend.position = "none"), gg_acoustic.overlap.group + theme_classic(base_size = bs) +
    theme(legend.position = "none"), gg_coeffs[[grep("count", names(gg_coeffs))]] +
    theme_classic(base_size = bs), gg_coeffs[[grep("area", names(gg_coeffs))]] +
    theme_classic(base_size = bs), gg_coeffs[[grep("overlap ~", names(gg_coeffs))]] +
    theme_classic(base_size = bs), gg_coeffs[[grep("group", names(gg_coeffs))]] +
    theme_classic(base_size = bs), nrow = 2, rel_heights = c(1.8,
    1))

# cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_70dpi.jpeg', width = 25,
# height = 9) cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_300dpi.jpeg', dpi = 300,
# width = 25, height = 9)

Takeaways

  • Higher vocal output in “high stress” birds compared to control

  • Higher acoustic overlap to themselves in week 1 for “high stress” birds compared to control

  • Decrease in overlap to themselves in week 1 across time

 


Session information

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] PhenotypeSpace_0.1.0 warbleR_1.1.26       NatureSounds_1.0.4  
##  [4] seewave_2.1.8        tuneR_1.3.3.1        ggridges_0.5.3      
##  [7] posterior_1.2.0      ggrepel_0.9.1        cowplot_1.1.1       
## [10] tidybayes_3.0.2      modelr_0.1.8         tidyr_1.2.0         
## [13] forcats_0.5.1        purrr_0.3.4          dplyr_1.0.8         
## [16] lme4_1.1-28          Matrix_1.3-4         brms_2.16.3         
## [19] Rcpp_1.0.8           GGally_2.1.2         sp_1.4-6            
## [22] MASS_7.3-54          formatR_1.11         knitr_1.37          
## [25] kableExtra_1.3.4     ggplot2_3.3.5        viridis_0.6.2       
## [28] viridisLite_0.4.0    pbapply_1.5-0        readxl_1.3.1        
## [31] lubridate_1.8.0      remotes_2.4.2       
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2            tidyselect_1.1.2      fftw_1.0-6.1         
##   [4] htmlwidgets_1.5.4     grid_4.1.0            munsell_0.5.0        
##   [7] codetools_0.2-18      DT_0.20               miniUI_0.1.1.1       
##  [10] withr_2.4.3           spatstat.random_2.1-0 Brobdingnag_1.2-7    
##  [13] colorspace_2.0-3      highr_0.9             uuid_1.1-0           
##  [16] rstudioapi_0.13       stats4_4.1.0          dtw_1.22-3           
##  [19] tensor_1.5            bayesplot_1.8.1       labeling_0.4.2       
##  [22] emmeans_1.8.1-1       rstan_2.21.3          polyclip_1.10-0      
##  [25] farver_2.1.0          bridgesampling_1.1-2  rprojroot_2.0.2      
##  [28] coda_0.19-4           vctrs_0.3.8           generics_0.1.2       
##  [31] TH.data_1.1-0         xfun_0.29             R6_2.5.1             
##  [34] markdown_1.1          gamm4_0.2-6           projpred_2.0.2       
##  [37] bitops_1.0-7          spatstat.utils_2.3-0  reshape_0.8.8        
##  [40] assertthat_0.2.1      promises_1.2.0.1      scales_1.1.1         
##  [43] multcomp_1.4-17       gtable_0.3.0          goftest_1.2-3        
##  [46] processx_3.5.2        xaringanExtra_0.7.0   sandwich_3.0-1       
##  [49] rlang_1.0.1           systemfonts_1.0.4     splines_4.1.0        
##  [52] spatstat.geom_2.3-2   broom_0.7.12          checkmate_2.0.0      
##  [55] inline_0.3.19         yaml_2.3.5            reshape2_1.4.4       
##  [58] abind_1.4-5           threejs_0.3.3         crosstalk_1.2.0      
##  [61] backports_1.4.1       httpuv_1.6.5          rsconnect_0.8.25     
##  [64] tensorA_0.36.2        tools_4.1.0           ellipsis_0.3.2       
##  [67] spatstat.core_2.4-0   raster_3.5-15         jquerylib_0.1.4      
##  [70] RColorBrewer_1.1-2    proxy_0.4-26          plyr_1.8.6           
##  [73] base64enc_0.1-3       RCurl_1.98-1.6        ps_1.6.0             
##  [76] prettyunits_1.1.1     rpart_4.1-15          deldir_1.0-6         
##  [79] zoo_1.8-9             cluster_2.1.2         magrittr_2.0.2       
##  [82] ggdist_3.1.0          colourpicker_1.1.1    mvtnorm_1.1-3        
##  [85] matrixStats_0.61.0    shinyjs_2.1.0         mime_0.12            
##  [88] evaluate_0.15         arrayhelpers_1.1-0    xtable_1.8-4         
##  [91] shinystan_2.5.0       gridExtra_2.3         rstantools_2.1.1     
##  [94] compiler_4.1.0        tibble_3.1.6          crayon_1.5.0         
##  [97] minqa_1.2.4           StanHeaders_2.21.0-7  htmltools_0.5.2      
## [100] mgcv_1.8-36           later_1.3.0           RcppParallel_5.1.5   
## [103] DBI_1.1.1             boot_1.3-28           permute_0.9-7        
## [106] cli_3.2.0             parallel_4.1.0        igraph_1.2.11        
## [109] pkgconfig_2.0.3       signal_0.7-7          terra_1.5-21         
## [112] spatstat.sparse_2.1-0 xml2_1.3.3            svUnit_1.0.6         
## [115] dygraphs_1.1.1.6      svglite_2.1.0         bslib_0.3.1          
## [118] webshot_0.5.2         estimability_1.4.1    rvest_1.0.2          
## [121] stringr_1.4.0         distributional_0.3.0  callr_3.7.0          
## [124] digest_0.6.29         vegan_2.5-7           spatstat.data_2.1-2  
## [127] rmarkdown_2.11        cellranger_1.1.0      shiny_1.7.1          
## [130] gtools_3.9.2          rjson_0.2.21          nloptr_2.0.0         
## [133] lifecycle_1.0.1       nlme_3.1-152          jsonlite_1.8.0       
## [136] fansi_1.0.2           pillar_1.7.0          lattice_0.20-44      
## [139] loo_2.4.1             fastmap_1.1.0         httr_1.4.2           
## [142] pkgbuild_1.3.1        survival_3.2-11       glue_1.6.1           
## [145] xts_0.12.1            shinythemes_1.2.0     stringi_1.7.6        
## [148] sass_0.4.0